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We present a novel joint time-dependent density-functional theory for the description of solute- 
solvent systems in time-dependent external potentials. Starting with the exact quantum-mechanical 
action functional for both electrons and nuclei, we systematically eliminate solvent degrees of free- 
dom and thus arrive at coarse-grained action functionals which retain the highly accurate ab initio 
description for the solute and are, in principle, exact. This procedure allows us to examine approx- 
imations underlying popular embedding theories for excited states. Finally, we introduce a novel 
approximate action functional for the solute-water system and compute the solvato-chromic shift 
of the lowest singlet excited state of formaldehyde in aqueous solution, which is in good agreement 
with experimental findings. 
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I. INTRODUCTION 

Electronic excited states are important in many ar- 
eas of physics, chemistry and biology. They are probed 
in spectroscopic experiments, such as absorption, re- 
flectivity or photoluminescence measurements. In addi- 
tion, they are important in many technical applications, 
such as photovoltaics 0, [|[ , laser technology [|| or light- 
emitting diodes[3, @. However, in most situations the 
excited system is not in isolation, but in contact with an 
environment. For example, dyes in Gratzel cells @, @] 
are in contact with an electrolyte and spectra of DNA 
molecules are typically obtained in an aqueous solution 
II- 

Various methods for the theoretical modelling of elec- 
tronic excited states have been developed. For extended 
systems, such as periodic solids or surfaces, excitation 
energies are typically extracted from the single-particle 
and two-particle Green's functions by solving the quasi- 
particle [10( and the Bethe-Salpeter equation [llj in the 
GW approximation. Quantum chemistry methods, such 
as configuration interaction 12] or coupled cluster the- 
ory [l3j], yield highly accurate excitation energies for 
atoms and small molecules. In contrast to the aforemen- 
tioned methods, which scale unfavorably with the system 
size, in recent years time-dependent density-functional 
theory [HGH emerged as an economical yet accurate the- 
ory for larger molecules and clusters d, [H, E3] • 

However, despite its good scaling properties, the ap- 
plication of time-dependent density-functional theory to 
electronic systems, which are not in isolation, but in 
contact with an environment, remains numerically chal- 
lenging. To capture solvent effects on excited states, 
a number of embedding approaches have been devel- 
oped (l8l - [2"2l | . These methods either model the solvent 
atomistically, for example by using the classical molec- 



ular dynamics technique [191 l2lj . or via a continuum 
approach [3 US H3- Due t° its simplicity, the latter ap- 
proach has enjoyed great popularity. In particular, many 
calculations employed time-dependent density-functional 
theory in conjunction with the "polarizable continuum 
model" [1, 0, E, HH , where the molecule is placed inside 
a cavity in a linear dielectric medium. The solute-solvent 
interactions are then separated into equilibrium and non- 
equilibrium contributions accounting for the fact that the 
electronic excitations on the molecule are screened by the 
high-frequency dielectric constant, which in many sys- 
tems is much smaller than the static dielectric constant. 



Despite the success of these continuum models, it is 
important to recall that their construction is purely phe- 
nomenological. To improve upon these theories and un- 
derstand their limitations, it is necessary to understand 
their origin from first principles. In this paper, we derive 
different levels of continuum embedding theory starting 
from the exact quantum-mechanical action functional for 
the full solute-solvent system (Sections [TT|. Next, as a 
demonstration that the approach can lead to practical 
calculations, we introduce a novel time-dependent con- 
tinuum model (Section IlIII) . which takes into account the 
deviation from bulk behavior of the solvent response in 
the first solvation shells and also retains the full frequency 
dependence of the dielectric response, and then apply this 
functional to the excitations of a formaldehyde molecule 
in aqueous solution ( Section ITVT) . Finally, in Section [V] 
we discuss our conclusions and describe possible future 
developments and applications. 
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II. JOINT TIME-DEPENDENT 
DENSITY-FUNCTIONAL THEORY 

In this section, we consider a physical system composed 
of an explicit subsystem (solute), in contact with an en- 
vironment (solvent). Both solute and solvent consist of 
electrons and nuclei of various species. Examples of such 
solute-solvent systems are a molecule dissolved in water 
or a defect in a host crystal. To investigate the behavior 
of a solute-solvent system in a time-dependent external 
potential, we employ time-dependent density- functional 
theory Here, we begin with the straightforward gen- 
eralization of the standard, fully causal expression for the 
action of a single-component system at zero temperature 
24] to multi-component systems at finite temperature. 
This gen eralization is similar to the theory of Li and 
Tong [23]. Those authors, however, worked only at zero 
temperature and employed the Frenkel-Dirac action func- 
tional, which violates the causality requirement of van 
Leeuwen [2~ij . Generalizing the approach of van Leeuwen 
to the present situation yields for the action functional S 
of the full solute-solvent system, 

S[n,{N a }} =A[n,{N a }}- j dr j drn{r, r)v(r, r) 

Jn Jc 

-Y, I dr I drN a (r,T)V a (r,T), (1) 
a Jn Jc 

where n(r) and N a (r) denote the density of electrons 
and the various nuclear species, while v(r) and V a (r) de- 
note the respective external potentials. Also, denotes 
an open volume as required when working in the grand 
canonical ensemble and r denotes the Keldysh time [26[ , 
which is defined on the contour C ranging from to oo 
just above the real time axis and then back from oo to 
just under the real time axis, and finally from to — ij3 
on the imaginary time axis, with /3 = l/(fcsT) being 
the inverse thermal energy. Finally, in ([T]), the intrinsic 
action A is the Legendre transform, with respect to the 
potentials v and V a [HI, of 



A[v,{V a }] = 
i log Tr < exp j3 



U(-i(3;0) 



(2) 



where /i c i and /i a are the chemical potentials of the 
electrons and nuclei, while 7V c i and N a denote the re- 
spective particle number operators. In the above, U = 
Tc exp(— i J cItH(t)) denotes the quantum-mechanical 
evolution operator [24| with Tc being the Keldysh time- 
ordering operator and H (r) denotes the standard many- 
body Hamiltonian for the electrons and nuclei of the 
solute-solvent system. 

We note that, relatively recently, Butriy et al. have 
also considered multicomponent time-dependent density- 
functional theory [27j |. but chose as the variational vari- 
ables the electron density in body-fixed coordinates and 



the diagonal of the nuclear N-body density matrix, 
whereas we here employ as the variables the much more 
tractable densities n(r,r) and N a (r,r). Butriy et al. 
employed their formalism to study correlated electron- 
nuclear excitations in isolated molecules. By contrast, we 
here are interested in electronic excitations of the system 
while treating the solute nuclei in the Born-Oppenheimer 
approximation, holding them fixed in place, so that they 
present a simple fixed external potential in which the elec- 
trons and solvent nuclei evolve. With the solute nuclear 
coordinates fixed, there is no need to work in body-fixed 
coordinates and a simple density description is sufficient 
for us to extract the density fluctuations of interest. 

We do, however, find it mathematically convenient as 
a matter of bookkeeping to treat the solute and environ- 
ment nuclear densities on an equal footing for as long as 
possible; therefore, we treat the solute nuclear densities 
as time-dependent in our derivation and only fix the loca- 
tions of the solute nuclei in the final step. Consequently, 
the index a in ([TJ and ([2]) ranges over all nuclear species 
in both the solvent and the solute. 

Because of the many environment degrees of freedom, 
finding the time-dependent densities which make the ac- 
tion in ([I} stationary is numerically challenging. More- 
over, the explicit details of the density fields describing 
the solvent are often irrelevant, because one is typically 
interested in properties of the solute. We therefore seek 
a fundamental description which treats the solute explic- 
itly and the solvent either at a simplified level or implic- 
itly. Petrosyan and coworkers j28|, [2!| have developed 
just such a rigorous "joint" density-functional theory for 
the static, equilibrium case. Specifically, Petrosyan et 
al. first minimize the full solute-solvent free-energy func- 
tional over the solvent electron density to arrive at a 
free-energy functional in terms of the solute electron and 
the solvent nuclear densities. The resulting theory treats 
both solute and solvent explicitly, but the solvent at a 
more coarse-grained level. Ultimately, for specific solute- 
solvent systems the coarse-grained free-energy functional 
is minimized over all solute electron and solvent nuclear 
densities to obtain the free energy of the overall system 
and its equilibrium properties. Petrosyan et al. devel- 
oped accurate and numerically tractable approximations 
to the coarse-grained functional and employed them to 
study surfaces and small molecules in aqueous solution 
with encouraging results (28l. [29j. 

To generalize the theory of Petrosyan and coworkers 
to the present non-equilibrium context, we split the total 
electron density in ([T]) into solute (n s ) and environment 
(n e ) contributions, n(r,r) = n s (r,r) +n e (r,r). Fun- 
damentally, a rigorous partitioning of electrons into so- 
lute and environment electrons is, of course, impossible 
because of the their quantum-mechanical indistinguisha- 
bility. Nonetheless, making S stationary with respect 
to all physically allowed environment electron densities 
and subsequently with respect to all physically allowed 
solute electron densities is guaranteed to recover the cor- 
rect total electron density. There are, of course, many 
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ways to express the total electron density as a sum of two 
subsystem densities. Thus, instead of a unique solution, 
there exists, in principle, a vast degenerate set of solu- 
tions in joint time-dependent density-functional theory 
consisting of all solute and environment electron densi- 
ties which sum up to the correct total electron density. 
In practice, however, we find that practical approxima- 
tions break this degeneracy and pick out a sensible par- 
ticular solution. This is reminiscent of the equilibrium 
case, where Petrosyan et al. 1291 observed that the use 
of molecular pseudopotentials [301 leads to sensible non- 



where stat„ e indicates that the expression in curly brack- 
ets is made stationary with respect to variations of n e 
and the superscript of A makes explicit the additional 
dependence of this functional on the external potentials. 
Note that we include the coupling term for the nuclear 
densities (%2N a V a ) in A: this partitioning is not func- 
tionally necessary, because the coupling term does not 
depend on n e and therefore maintains its simple form. 
This choice, however, ensures that A describes a neu- 
tral system: Maintaining charge neutrality is important 



where the first four terms describe the solute: The first 
term, A e i^ s , denotes the intrinsic action of the solute elec- 
trons and is typically [24| written as 

A c i, s = A KS -A h - A XC , (5) 

where Aks denotes the action of non-interacting elec- 
trons, Ah = 1/2 J dr J dr' J dTn s (r,T)n s (r',T)/\r — r'\ 
is the Hartree contribution and Axe the exchange- 
correlation term. The second and third terms in (QJ, 
respectively, are the intrinsic action of the solute nuclei 
(with densities N s , a ) and their coupling to the exter- 
nal potentials. The fourth term captures the interac- 
tion between solute electrons and solute nuclei. In our 



degenerate solutions. 



A. Explicit solvent functionals 

Making the action stationary with respect to the en- 
vironment electron density while holding the solute elec- 
tron and environment nuclear densities fixed, we obtain 
the coarse-grained explicit-solvent functional 5 CX , 



,r)+Y,N a (r,T)V a (r,T) 




(3) 



both formally to ensure the existence of the thermody- 
namic limit [31| and also practically to mitigate the need 
to capture long-range couplings within an approximate 
functional. We stress that our partitioning does not fun- 
damentally complicate the functional dependence of the 
new functional A because the coupling to V a will retain 
its simple form in terms of N a . 

To find practical approximations, we partition A into 
various physically meaningful contributions according to 



(4) 

I 

actual calculations, we hold the solute nuclei fixed in 
space, finding A nuc/cl s = / dr J dTn s (r, r)w st (r), with 
v a t being the static potential created by the solute nu- 
clei. Under these conditions, the second and third term 
in Q become independent of all time-dependent degrees 
of freedom (though, not of time-dependent potentials) 
and can be dropped for the purpose of the variational 
calculations. However, for later convenience we retain 
Amc, s = / dTj2i<j ZiZj/\Ri - Rj\ with Ri and Zj 
denoting the positions and charges of the solute nuclei. 

The fifth term in describes the isolated neutral envi- 
ronment in terms of its nuclear densities, N e<a . Accurate 
approximations for this action functional are less well 



S ex [n s , {N a }} =stat„ e < A[n s + n e , {N a }} 



dr / dr 

Q JC 



n e (r,T)v(r 



dr / dTn s (r,T)v(r,T) 
n Jc 

A Mv ^[n s ,{N a }}- f dr f dTn s (r,r)v(r,T), 



r 



i 

A =A chs [n s ] + A nuc , s [{N s , a }] - A^^ s [{N Sta }} - A nuc/el)S [n s , {7V S , Q }] 

+ A e [{N e , a }] - A^-^JV^}] - A^K(v;,.Mv. 1< ,}) [naj {Ns>a}j {Ne>a}l 
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known than for electrons. However, there has been much 
progress recently in the construction of such action func- 
tionals for classical liquids [H, , which constitute an 
important and technologically relevant class of solvents. 

Finally, having identified all interactions between 
charged species, we expect the remaining two contri- 
butions in (0| to be relatively small. The sixth term, 
A^4 C xt,c, describes the interaction between the neutral 
solvent and the external potentials. The final, seventh 
term, AA s . e , maintains the full functional dependence 
and captures by definition all remaining interactions. We 
find that the coupling between the neutral solute and the 
neutral solvent constitutes the most important contribu- 
tion to this term. 

To find approximations to A_4 0X t lG , we observe that the 
best form for this term depends on the physical system 
under consideration, because the solvent electrons screen 
the bare nuclear charges in qualitatively different ways 
depending on the physical nature of the solvent. If the 
solvent consists of ions of charge Z a , the corresponding 



action would be ^ Q / dr J dr[N etCt V a + (Z a — Z a )N e ^ a v] 
with Z a being the true charge of the nucleus. If, how- 
ever, the solvent is composed of neutral polar molecules, 
where each "effective" nucleus carries a partial charge 
q a and all partial charges in a molecule add up to 
zero, one should replace Z a in the above expression by 
q a . Finally, if the solvent consists of apolar molecules 
or neutral atoms, we can approximate the coupling by 
J2 a J dr J drir a N e ^ a \Vv\ 2 with n a being the polarizabil- 
ity of the "effective" nucleus a. 

The final term, AA s , e , has the full functional depen- 
dence and thus can capture all remaining interactions. 
In typical density-functional theory fashion, because we 
have separated out by various approximations all other 
possible interactions and ensured that this term repre- 
sents a charge-neutral interaction, we expect this term 
to be relatively small with mild functional dependencies 
and thus amenable to simple approximations. Accord- 
ingly, we expand AA Sj e as a Taylor series in the various 
densities, keeping only the lowest-order coupling terms, 



AAs.e = [ dr [ dr' [ dr [ dr' V N e , a (r', r') 
Jn Jn Jc Jc 



u a (r, r', r, r')n s (r, r) + ^ w ap{r, r', r, t')N s ^ (r', r') 

P 



(6) 



where w a (r,r',T,T') = 5 2 AAs.e/5n s (r,T)5N e . a (r',T') 
and w a p(r,r',T,T r ) = 5 2 AA s , e /SN e , a (r,T)SN s ^(r',T') 
denote effective time-dependent interaction potentials 
between solvent nuclei and solute electrons or solute 
nuclei. Note that, in principle, the Taylor series con- 
tains various other coupling terms: for example, a term 
quadratic in the solute electron density can occur. How- 
ever, such a term only renormalizes the Hartree contri- 
bution in .Aci.s and is therefore neglected in ([6]). 

Tablc|I]summarizes all of the above considerations, list- 
ing all of the various contributions to S C x and the terms 
which capture them. 



B. Implicit solvent functionals 

1. General considerations 

Rather than follow the above route of dealing explic- 
itly with the solvent nuclei, for this initial work, we take 
a simpler tack that allows us to make contact with stan- 
dard continuum solvent models. For this purpose, we 
eliminate the environment nuclei from S ex and introduce 
a new action functional <Si m which depends on the solute 



densities only and treats the solvent implicitly as follows, 
S im [n s ,{iV s , a }] = stat {Neia} A^ v ^[n s ,{N s , a },{N e , a }} 

— dr drn s (r,T)v(r,T) 
Jn Jc 

= 0(«'{ v «»[n„{iV. |a }]- f dr [ dTn s (r,T)v(r,T). 

Jn Jc 

(7) 

Again, we partition Q into meaningful contributions ac- 
cording to 

G = A el , s [n s ] + A nuc , s [{N s , a }} - A^/HiiNs,.}] 

- A nuc/eU [n s ,{N s , a }] - AG^ v ^ln s ,{N s , a }} (8) 

with AQ = -stat {A r e , a }[A e - A_4 cxti0 - AA s , e ] and _4 c i, s 
is given by ([5]). Note that AQ depends on the solute 
densities of electrons and nuclei, but also on the time- 
dependent external potential. To understand the conse- 
quences of the additional functional dependency, we now 
investigate the linear response behavior of S{ m in greater 
detail (assuming fixed solute nuclei). 

The time-dependent solute electron density corre- 
sponding to v makes Si m stationary, SSi m = 0, which 
implies 8Q/8n s = v. Using ([H]) we thus have (assuming 
fixed solute nuclei) 

v(r,r) =v KS (r,T) - v H (r,r) - v X c{r,r) 

-« st (r)-^-»(r,r), (9) 
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TABLE I: The Table shows the various contributions to 
Sox due to interactions between solute nuclei (N B ), solute 
electrons (n s ), solvent nuclei (N e ), solvent electrons (n e ) 
and the external potentials (v and V). Note that in S ex the 
solvent electrons are treated implicitly. 



find 



contribution 


contained in 


Ns/Ns 




N s /n s 


■^nuc / cl,s 


N s /N e 


AAs,e 


N s /n e 


AAs,e 


N s /V 


A 

^nuc/cxtjS 


n s /n s 




n s /N e 


AAs,e 


n s /n e 


AAs,e 


n a /v 


last term in ([J 


N e /N e 




N e /n e 


A e 


N e /V 




n e /n e 


A e 


Tle/v 





8vKs{ r , T ) — Sv(r, t) + / dr 



,6n 8 (r',T) 



+ [ dr' [ dT'f X c(r,r',T,T')5n s (r',T') 
Jn Jc 

+ f dr' f dT'ff(r,r',T 7 T')5n s (r',T') 
Jn Jc 

+ f dr' f dT'f*(r,r',T,T')5v(r',T') 
Jn Jc 

+ E / dr ' I dT'fZ a (r,r',T,T')5V e , a (r',T'), 
a Jn Jc 

(11) 

where ff = Sv e /Sn s , ff = Sv e /Sv and / e c Q = Sv e /SV e , a 
denote additional contributions to Svks caused by the 
environment. In an actual experiment, where the whole 
solute-solvent system is probed, for example by an elec- 
tromagnetic wave, we expect SV e , a to be related to 6v. 
In this case, we can express the last term in (fTTj) as fj?8v 

with/,? = E a /, c A./k 

The change in the Kohn-Sham potential is related to 
(5n s via 



where vks — SAKs/Sn s denotes the Kohn-Sham po- 
tential, v H = SA H /6n s , v X c = SA X c/8n s and v st 
is the static potential due to the solute nuclei. Also, 
v e = (5 AC? / 8n s denotes the additional potential due to the 
presence of the environment. Note that AQ and there- 
fore also v e depend both on the solute density n s and the 
external potential separately. 

A small change Sv in the external potential causes a 
change Sn s in the solute electron density. In the linear 
response regime, these quantities are related via the re- 
sponse function \ 



6n s (r,r)= [ dr' j dr'x(r,r' ,T,T')6v(r' ,t'). (10) 
Jn Jc 



To compute x, which is the observable in spectroscopic 
experiments on the solute-solvent system, we first deter- 
mine the change in the Kohn-Sham potential Svks cor- 
responding to 8v. (Strictly speaking, only the total re- 
sponse function Xtot — Sn/5v = Sn s /6v + 5n e /5v is mea- 
sured. However, for solute-solvent systems where the re- 
sponse of the solute occurs in a different frequency range 
than the response of the solvent one can determine ex- 
perimentally the solute response function x- This is the 
case for the lowest singlet excitation of formaldehyde in 
water, which we study in the Section ITVl ) Using ©, we 



5n s (r,r)= / dr' dr'xKs(r,r',r,T')SvKs(r',T'), 
Jn Jc 

where xks denotes the response function of non- 
interacting electrons. Combining (|12l) . (fTT| and (fTU|) and 
adopting a matrix formulation for the space and Keldysh- 
time variables then yields 

= [1 + fe + fe}' 1 {XKS -\K+ IXC + ft] } , 

(13) 

where K denotes the matrix corresponding to the 
Coulomb interaction K (r, r', t, t') = S(t, T')/\r— r'\ with 
8{t, t') denoting the Delta- function on the Keldysh con- 
tour. 

Compared to the familiar equation for x without sol- 
vent, given by x 1 = Xks ~ [K + fxc], we find that 
(|13[) contains three extra terms due to the presence of 
the solvent: describes the change of the solvent po- 
tential due to a change in the solute electron density, 
while ff and /? describe changes induced by a vari- 
ation in the external potential. Without justification, 
"polarizable continuum model" approaches typically [Tij 
approximate the potential due to the solvent as a func- 
tional of the solute density only, which means they only 
include ff and neglect ff and /P . This insight into the 
assumptions underlying popular embedding approaches 
underscores the value of following the density-functional 
approach rigorously so as to identify all potentially rele- 
vant functional dependencies. 
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2. Practical approximations 

To develop practical approximations, we separate AQ 
into a contribution AQ se [n SJ {iV S Q }], which describes the 



where W(r, r', r, r'), W a (r, r' , r, r') and W a p(r, r', r, r') 
denote effective interaction potentials between the vari- 
ous solute particles. In the next section, we approxi- 
mate these interaction potentials by screened Coulomb 
interactions, which results in a simple, yet accurate joint 
density-functional theory for solute-water systems. 

Approximating AQ cxt is more difficult: a possible route 
to finding explicit functionals is to express the environ- 
ment nuclear densities in terms of solute densities ac- 
cording to N eyCC (r) = g a [n s , {N Sy p}](r) and insert this 
relation into the various forms for AA cx t,c discussed in 
the last section. We expect that g a has a similar form 
as the dielectric function employed in the next section 
[Equation (|18[)]. where we employ a local ansatz to de- 
scribe the crossover from bulk screening to vacuum. 

However, for our first implementation of joint time- 
dependent density-functional theory for solute- water sys- 
tems presented in the next section we neglect AQ cxt ■ We 
expect, however, that the additional solvent response due 
to this term can be included by "renormalizing" the di- 
electric function describing the environment (see Section 
IIIip . Future work should explore the consequences and 
importance of this term. 

III. AN IMPLICIT ACTION FUNCTIONAL FOR 
THE WATER-SOLUTE SYSTEM 

To allow us to explore and test the potential of the 
above ideas in an actual application, in this section, we 



with p s (r,r) — —n s (r,T) + J2i Zi5(r — Ri) denoting the 
solute charge density and K = K + W, where K is the 
bare Coulomb interaction defined as above. Also, AV ps 



interaction among solute particles mediated by the envi- 
ronment, and a remainder AQ^^ Ve ' [n s , {A SjQ }]. Tay- 
lor expanding AQ s e yields 



(14) 



introduce a relatively simple, approximate joint time- 
dependent density-functional for the solute- water system. 
In particular, we assume that all solvent effects can be 
described via a position and frequency-dependent local 
dielectric function, which depends on the electronic struc- 
ture of the solute. Inclusion of the spatial dependence of 
screening effects is crucial, because the dielectric response 
of water in the first solvation shells differs notably from 
the bulk response. Also, in contrast to "polarizable con- 
tinuum model" approaches [l8l . l34j |. where a particular 
value for the high-frequency dielectric constant is cho- 
sen, we employ the full frequency-dependent dielectric 
function. 



Specifically, the assumption of dielectric screening im- 
plies that all effective interactions introduced in (|14l) are 
proportional to the screened interaction W between two 
unit charges and only rescaled by the charges of the inter- 
acting species. In particular, we approximate W = W, 
W a = — Z a W and W a p = Z a ZpW. The resulting action 
functional for the solute- water system is then given by 



(15) 

I 

reflects the fact that, in practical calculations, we employ 
the pseudopotential approximation (35j |. in which the nu- 
clei are replaced with ionic cores of charge Zj, whose 



Ag s , e = 



dr / dr' dr dr' 

n Jn Jc Jc 



-n s (r, r)W(r, r' , r, T r )n s (r', r') 



J 



Q = A KS - A XC - AV P s - \ [ dr [ dr' [ dr [ dr' p s (r,T)K(r,r' ,r,r)p s (r' ,r), 

^ Jn Jn Jc Jc 
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potentials at large distances (when not screened by the 
environment) go as Zj/\r — Rj\ but which differ from 
this by a localized function AV ps (i — Rj) within a small 
"core radius" that represents a distance much smaller 
than where we would expect screening from the environ- 
ment to occur. Within our framework, the long-range 
parts enter through the solution of (|16p and thus are 
properly screened, and the short-range parts contained 
in AVp S enter directly as they require no such screening. 

The screened potential corresponding to a physical 
charge density p s (r,t), which is equal on both vertical 
branches of the Keldysh contour, is given by 4> s = K R p s , 
where K R denotes the retarded interaction [24]]. In ac- 
tual calculations, we obtain (f> s by solving the screened 
Poisson equation 



Clausius-Mossotti form 



V • e(r,cj)V0(r,w) = — 47r ( o s (r,w). 



(16) 



All information about the environment is contained in 
the dielectric function e(r, uS). In principle, both the ionic 
and the electronic degrees of freedom of the solvent con- 
tribute to the dielectric response. We demonstrate below 
that for the frequencies of interest, we can safely ignore 
the motion of the ions and only deal with the electronic 
response corresponding to a fixed nuclear solvent density. 
We make the natural assumptions that the system is in 
equilibrium before the excitation and that the equilib- 
rium nuclear solvent density is determined locally by the 
equilibrium solute electron density rioir). This suggests 
the following local ansatz for the dielectric function, 



e(r,ui) = e(no(r),w). 



(17) 



This ansatz is physically reasonable in that it interpo- 
lates smoothly between the dielectric response of vacuum 
and the bulk liquid and thus avoids the need to specify 
a cavity shape. If we further assume that the frequency 
dependence of e(no(T"),w) enters only through the fre- 
quency dependence of the bulk dielectric function €& (<*;), 
we can generalize the form employed by Petrosyan and 
co-workers I28| to 



e{r,u}) = l-\ erfc 



log(n (r)/n c 
V2a 



(18) 



where the parameters n c and a determine the location 
and width, respectively, of the crossover from the vac- 
uum to the bulk liquid dielectric response. Petrosyan et 
al. [28[ determined the numerical values n c = 4.73 x 
10 -3 A -3 and a = 0.6 for these parameters by fitting 
solvation energies of small molecules obtained by their 
equilibrium joint density-functional theory to experimen- 
tal data. We choose to work with these values as well. 

To complete the theory, we need an expression for the 
frequency-dependent bulk dielectric constant At 
frequencies corresponding to electronic excitations, we 
may ignore the complicated low-frequency dielectric re- 
sponse of water and employ a model which describes 
the high-frequency range reliably. For this, we use the 



e b (ui)-l 47r 

^HT2 = T^H, 



(19) 



where n b denotes the bulk molecular particle density of 
water and a(ui) = J2j Fj/(Ej~ uj2 ) denotes the mean po- 
larizability of an isolated water molecule, with Fj and Ej 
being the oscillator strength and excitation frequencies, 
respectively, for excited state j. In the next section, we 
compute a(uj) using time-dependent density-functional 
theory and demonstrate that flT9"|) . which neglects the 
contribution from the permanent dipole moments, indeed 
reliably describes the bulk screening response of water to 
low-lying electronic excitations. To obtain excitation en- 
ergies of the solute, we analyze the linear response of (fT5|) 
resulting in 



X" 1 = Xks ~ [K + fxc], 



(20) 



which lacks the subtleties appearing in Section III Bl be- 
cause the present model lacks any explicit environment 
dependence on the external potential. Our final work- 
ing equation is obtained by expressing (|20l) in transition- 
space notation[TH. where the fused index k — (k,j) de- 
notes a transition between two equilibrium Kohn-Sham 
orbitals ipj(r) and ib^r). We arrive at a self-consistent 
eigenvalue problem [151 ] for excitation energies Ej of the 
solute, 



(21) 

where Ae K = ej — tk with ek denoting the equilibrium 
orbital energies and the eigenvector C« determines the 



oscillator strength of the transition If 
matrix is given by 



The coupling 



dr / dr'$*(r) K R (r,r',cj) + f£ c (r,r',uj) $„(r') 
Q ill L J 

(22) 



with 3> re (r) = 4'k( r )' t l J j( r ) an d fxc denotes the re- 
tarded exchange-correlation kernel. Note that even for a 
frequency-independent exchange-correlation kernel, the 
solvent response makes M frequency- dependent. Equa- 
tion (l2"Tj) is solved iteratively: setting K R (Ej) = K 

yields an initial estimate E^ for the excitation energy. 

Next, we solve (|2"Tj) using K R {E^) and iterate until self- 
consistency is achieved. 



IV. APPLICATION TO FORMALDEHYDE IN 
AQUEOUS SOLUTION 

As a test case, we study the lowest singlet excited 
state of a formaldehyde molecule in aqueous solution. 
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A number of theoretical approaches have been applied 
to stud y so lvato-chromic shifts of formaldehyde in water 
flU [l9l. l34l. I37M40I ] . However, the agreement with experi- 
mental findings has generally been unsatisfactory. 

In this section, we first compute the mean polarizabil- 
ity of an isolated water molecule from time-dependent 
density-functional theory and obtain the frequency- 
dependent bulk dielectric function of liquid water using 
the Clausius-Mossotti equation. Next, we explore the ex- 
citations of formaldehyde in the gas phase and in solution 
using the joint time-dependent density-functional theory 
described in the last section. 

All calculations are carried out in a plane wave basis 
with a cutoff of 40 hartree. We use Kleinman-Bylander 
pseudopotentials [4l[ and a cubic supercell of length 
20 bohr. For the ground state calculations we employ 
the local density approximation [42], [43[ and for the ex- 
citations the adiabatic local density approximation [HI ]. 



TABLE II: Comparison of the lowest singlet excitation ener- 
gies of an isolated water molecule with previous theoretical 
work by Bernasconi [44| and experiment [44[- All results are 
given in eV. 



This work 


Ref.[44] 


Expt. [44] 


6.47 


6.39 


7.4 


7.74 


7.78 


9.1 


8.01 


8.05 


9.7 



TABLE III: Comparison of our results for the static mean 
polarizability of an isolated water molecule and for the 
optical dielectric constant e opt = et(iJ = 1 eV) of liquid 
water with previous theoretical work by Hu et al. [4(| and 
experiment [46H48T ]. 



A. Dielectric function of liquid water 

To compute excitation energies of a formaldehyde 
molecule in aqueous solution, we need the frequency- 
dependent dielectric response of water. According to 
(fin]), this requires the mean polarizability of an isolated 
water molecule. We first carry out ground-state calcula- 
tions and fully relax the electronic and ionic structure. 
Then, we employ time-dependent density-functional the- 
ory to obtain excitation energies and oscillator strengths 
using all 4 occupied Kohn-Sham orbitals, plus an addi- 
tional 220 unoccupied orbitals. Table |TT] shows our re- 
sults for the three lowest singlet excitation energies of an 
isolated water molecule, and compares them to previous 
theoretical work [lij and also to experiment [44[. The 
discrepancy between theory and experiment is around 
1 eV or larger for all excited states. The poor perfor- 
mance of time-dependent density-functional theory for 
the water molecule can be traced to the Rydberg char- 
acter of the excitations, which cannot be described in 
the adiabatic local density approximation due to the in- 
correct asymptotic behavior of the exchange-correlation 
potential at large distances 0, |H[ . 

Despite these problems, the adiabatic local density ap- 
proximation gives good results for the static polarizabil- 
ity and for the low-frequency dielectric constant of liquid 
water: Table [TTT] compares our results for these quantities 
with previous calculations [46[ and also with experiment 
[46M481 | . We observe that the Clausius-Mossotti formula 
(|19[) describes the dielectric response of liquid water very 
well in the frequency range corresponding to low-lying 
electronic excitations. 



B. Formaldehyde in the gas phase 

Next, we explore the excitations of formaldehyde in the 
gas phase. Table IIVI compares our results for the three 



Units This work Ref.[46] Expt. [46-48] 



5(w = 0) bohr 3 


10.50 


10.52 


9.6-9.9 


£opt 


1.83 




1.78 



lowest singlet excited states with a previous calculation 
by Bauernschmitt and Ahlrich [49( , who also employ the 
adiabatic local density approximation, and also with ex- 
perimental findings [18l . |49( . We observe that our excita- 
tion energy for the lowest state is relatively close to the 
experimental value, while the higher states deviate more 
than 1 eV from experiment. Again, the relatively large 
deviation for the higher excited states is due to the incor- 
rect long-distance behavior of the exchange-correlation 
potential in the adiabatic local density approximation. 

Comparison of Tables IIVI and [TT] shows that the lowest 
excitation energy of formaldehyde is several eV smaller 
than the corresponding value for water. Therefore, in 
our joint time-dependent density-functional calculations 
we only evaluate the dielectric function at frequencies 
smaller than its first pole. In this region, e&(w) is close 
to unity and can be approximated by a constant. This 
is a common approximation in "polarizable continuum 



TABLE IV: Comparison of the three lowest singlet excitation 
energies of an isolated formaldehyde molecule with previous 
theoretical work by Bauernschmitt et al. [4^] and experiment 
[T^.|49|. All energies are given in eV. 



This work 


Ref.[49j 


Expt. [18, 49] 


3.66 


3.64 


3.8-4.2 


5.68 


5.93 


7.13 


6.78 


6.79 


8.14 
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TABLE V: Comparison of our results for the equilibrium 
dipole moment of formaldehyde in vacuum (p V ac) and aqueous 
solution (paoiv) with experiment [H^] and previous theoretical 
work [331 . Dipole moments are given in bohr. 



blue-shift, we express the excitation energy as 





This work 


Ref. [38] 


Expt. [50] 


Pvac 


0.91 


0.90 


0.91 


Psolv 


1.32 


1.32 





TABLE VI: Lowest singlet excitation energy of formaldehyde 
in vacuum (E vac ) and aqueous solution (E so \v) obtained 
from joint time-dependent density-functional theory. The 
resulting solvato-chromic shift (i? so i v — _E V ac) is shown in the 
third column and compared to experiment [l!|. All results 
are given in eV. 





-f^solv 


Shift 


Expt. [19] 


3.66 


3.83 


0.17 


0.21 



model" approaches, which for the non-equilibrium re- 
sponse employ a frequency-independent dielectric func- 
tion derived from the index of refraction of water fl8| . 

However, this approximation breaks down for solutes 
with higher-energy excited states that are comparable or 
larger than the lowest excited state of water. In this case, 
the full frequency dependence of the dielectric response 
must be retained and a self-consistent solution of (j2"Tj) is 
necessary. 



C. Formaldehyde in aqueous solution 

We now apply joint time-dependent density-functional 
theory to calculate excitations of a solvated formalde- 
hyde molecule. We use the static joint density-functional 
theory of Petrosyan et al. [28j to determine the equilib- 
rium electronic structure neglecting ionic relaxations in- 
duced by the aqueous environment, which as was shown 
by Kongsted et al. only lead to shifts in the excitation 
energies of about 0.01 eV [19(. Table IVl shows that the 
equilibrium dipole moment obtained in our calculation is 
in excellent agreement with previous theoretical work [38$ . 

We then solve the linear response equation (f2"Tj) 
of joint time-dependent density-functional theory self- 
consistcntly, as described in the last section. The lowest 
excitation energy is converged to within 0.01 eV after two 
iterations. Table IVT1 summarizes our results for the exci- 
tation energies in vacuum and solution and the resulting 
solvato-chromic shift, which is in good agreement with 
the experimental value [l9|. 

To physically understand the observed solvato-chromic 



-^vac/sol V Aty 



c/solv Tvac/sol v > 



(23) 



where Ae vac / so i v denotes orbital energy difference in vac- 
uum or solution and 7 va c/soiv denotes the correction from 
(joint) time-dependent density-functional theory. For 
gas-phase formaldehyde, Ae vac = 3.34 eV already gives 
a reasonable approximation to E vajC = 3.66 eV with the 
correction 7 vac = 0.32 eV being relatively small. In so- 
lution, we find Ae so i v = 3.56 eV and 7 so iv = 0.27 eV. 
We conclude that the solvato-chromic shift is determined 
mostly by the change of the orbital energy differences, 
Ae vac — Ae so i v = 0.22 eV which is quite close to the to- 
tal solvato-chromic shift of 0.17 eV, while the correction 
term changes relatively little (only 0.05 eV). 

The change of Ae upon solvation is caused by the dif- 
ferent coupling of ground and excited states to the aque- 
ous environment: the ground state has a large dipole 
moment (j3gg C = 0.91 bohr) and couples strongly to the 
aqueous environment leading to a large negative solva- 
tion energy, while the dipole moment of the excited state 
(Pex° = 0-14 bohr) and the resulting solvation energy are 
much smaller. 

Table IVlTl compares our result for the solvato-chromic 
shift of the lowest singlet excitation with previous the- 
oretical work [H, [34], l37M40l ] and also with experiment 
19]. In contrast to our calculation, which gives very good 
agreement with the experimental data, calculations em- 
ploying the "polarizable continuum model" to describe 
the solvent underestimate the shift [HI, HH, while ap- 
proaches, which treat the solvent atomistically within a 
supermolecular approach [H, |39[, typically overestimate 
the shift. We point out that both methods suffer from 
weaknesses which are absent in our approach. In partic- 
ular, the results of the "polarizable continuum model" 
approach depend sensitively on the chosen cavity size 



and shape [5l| . This indicates that a more realistic de- 
scription of the solvent response in the vicinity of the so- 
lute is of great importance. Atomistic solvent models, on 
the other hand, offer a reliable description of the solvent 
structure close to the solute, but to compute converged 
thermodynamic averages the sampling of many solvent 
configurations is required. In addition, supermolecular 
approaches which model the solvent by a finite cluster 
surrounding the solute do not capture the long-range di- 
electric response of the solvent. 

Our approach includes both long-range screening ef- 
fects and a reliable description of the solvent response 
close to the solute. Similarly, Naka et al. [40| . who em- 
ploy the reference interaction-site model to describe the 
solute, and Kongsted et al. [371 ] . who combine an atom- 
istic treatment of the first solvation shells with the "po- 
larizable continuum model," also obtain solvato-chromic 
shifts in good agreement with experiment. However, un- 
like our action functional, these models are not derived 
from first principles: instead they start out with a par- 
tioning of the action and then approximate each con- 
tribution typically with a different level of theory mak- 
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TABLE VII: Comparison of our joint time-dependent 
density- functional theory results, previous theoretical work 
[Isl . Hj, l37l - |4Ci| and experiment pjj for the solvato-chromic 
shift of the lowest singlet excited state of formaldehyde in 
aqueous solution. The second column lists the solvent model: 
joint time-dependent density- functional theory (JTDDFT), 
the "polarizable continuum model" (PCM), the supermolec- 
ular approach (SM) or the reference interaction-site model 
(RISM). All results are given in eV. 





Method 


Shift 


Ref. [18] 


PCM 


0.12 


Ref. [34] 


PCM 


0.12 


Ref. [38] 


SM 


0.33 


Ref. [39] 


SM 


0.39 


Ref. [40J 


RISM 


0.25 


Ref. [37] 


SM+PCM 


0.23 


This work 


JTDDFT 


0.17 


Expt. [191 




0.21 



ing it difficult to judge the limits of their applicability 
a priori and to systematically improve upon them. For 
example, Naka et al. combine a CASSCF treatment of 
the solute electrons with an electrostatic coupling scheme 
between solute and solvent and the reference-interaction 
site model for the solvent structure [io1 |. 

V. SUMMARY AND CONCLUSIONS 

In sum, we describe the construction of a joint time- 
dependent density-functional theory for the modelling of 
solute-solvent systems. We derive coarse-grained action 
functionals by eliminating environment degrees of free- 
dom. This procedure enables us to examine the under- 
lying assumptions and uncover previously ignored func- 



tional dependencies in popular approaches such as the 
time-dependent "polarizable continuum model" and to 
explore their domains of validity. In particular, we find 
additional contributions to the action functional which 
are typically neglected in standard approaches. Also, in 
order to replace the full frequency dependent solvent re- 
sponse by a high-frequency dielectric constant, as is often 
done in standard approaches, the excitation energy of the 
solute has to be far from the poles of the solvent dielec- 
tric function. Otherwise, a self-consistent solution of the 
linear response equation is necessary. 

We also introduce an explicit, approximate action func- 
tional for the modelling of electronic systems in aque- 
ous solution. Application of this functional to solvated 
formaldehyde leads to good agreement with experiment 
for the solvato-chromic shift of the lowest singlet excited 
state. The novel implicit functional we introduced can 
now be applied to more complicated systems, such as 
dyes in Gratzel cells 0] or solvated DNA molecules @, Q . 

The framework is now in place for future work to de- 
velop approximate forms for the explicit solvent function- 
als by making use of existing forms of time-dependent 
functionals for classical liquids 33] and to generalize 
time-independent functionals for molecular liquids devel- 
oped by us [H, IE3 to include time-dependent response 
of the solvent electrons. The resulting theory would then 
allow us to describe the nuclear dynamics of the envi- 
ronment during the excitation and to compute solvation 
relaxation functions which have been measured experi- 
mentally [54 l55l|. 
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